arXiv: 1506.0591 Ivl [q-fin.PR] 19 Jun 2015 


Seasonal Stochastic Volatility and Correlation together with the 
Samnelson Effect in Commodity Fntnres Markets * 


Lorenz Schneider Bertrand Tavin ^ 


Jnne 22, 2015 


Abstract 

We introduce a multi-factor stochastic volatility model based on the CIR/Heston volatility process 
that incorporates seasonality and the Samnelson effect. First, we give conditions on the seasonal 
term under which the corresponding volatility factor is well-defined. These conditions appear to be 
rather mild. Second, we calculate the joint characteristic function of two futures prices for different 
maturities in the proposed model. This characteristic function is analytic. Finally, we provide 
numerical illustrations in terms of implied volatility and correlation produced by the proposed model 
with five different specifications of the seasonality pattern. The model is found to be able to produce 
volatility smiles at the same time as a volatility term-structure that exhibits the Samnelson effect 
with a seasonal component. Correlation, instantaneous or implied from calendar spread option prices 
via a Gaussian copula, is also found to be seasonal. 
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1 Introduction 


Seasonality is a well-known empirical feature of several commodities markets. In the energy sector, 
among fossil fuels, natural gas futures curves, and among refined products, gasoline, heating oil and fuel 
oil futures curves all typically display seasonality. In the agricultural sector, almost all futures curves 
show seasonality due to harvest times and the seasons of the year. 

It is important to distinguish from the outset between two types of seasonality: seasonality of futures 
prices and seasonality of volatility of futures prices. 

Regarding seasonality of prices, consider agricultural commodities such as corn, soybeans and wheat. 
These tend to be in high supply after the harvest in summer, and in low supply in the months preceding 
the harvest. This typically leads to relatively low prices of futures contracts with delivery months in the 
summer or early fall, and high futures prices of contracts with delivery months in late winter or spring. 

*We would like to thank Ennio Fedrizzi, Frangois Le Grand and Cassio Neri for helpful and stimulating comments, 
discussions and suggestions. 

^Center for Financial Risks Analysis (CEFRA), EMLYON Business School, schneiderOem-lyon.com. 

^Center for Financial Risks Analysis (CEFRA), EMLYON Business School, tavinOem-lyon.com. 
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Therefore, when the prices of these contracts are plotted as a function of their maturity, they tend to 
rise and fall with the maturity in some seasonal way. In other words, the futures curve shows seasonality. 
However, the price of an individual futures contract with a given maturity should not rise and fall over 
time in any kind of seasonal way: indeed, such a behaviour would lead to easy arbitrage opportunities. 


Regarding seasonality of volatility, the situation is different in the sense that now an individual futures 
contract, with fixed maturity, tends to go through phases of relatively high or low volatility according to 
a seasonal pattern. To take again the example of agricultural commodities, the weather in the months 
leading up to the harvest has a direct impact on its quality and quantity, and futures prices can fluctuate 
strongly as forecasts for the new crop change. In contrast to this, weather patterns in winter tend to be 
of minor consequence for the harvest, and futures prices tend to fluctuate less strongly. 


It follows from these empirical observations that for commodity models, seasonality is usually only an 
issue for the volatility, but not for the futures price itself. Mathematically, individual futures prices are 
model l ed as martingales, and martingales do not have a tendency to rise or fall in a pre-determined way. 
Clark ( 2 OI 4 II gives a general discussion and numerous examples of seasonality in various commodities 


markets. 


Traditionally, there are two approaches to modelling the prices of futures contacts: futures-based 
models and spot-based models. An advantage of futures-based models models is that since the futures 
price curve is an input of the model, any arbitrage-free shape of the initial futures curve can be accommo¬ 
dated, including any type of seasonality. In contrast, a first step for spot-based models is to make them 
fit the initial futures curve, which uses up model parameters and doesn’t necessarily lead to satisfactory 
results. 


[s^ renseni (|2002h studies the modelling of seasonality in corn, soybean and wheat futures markets. 
Analysis of a large data set of CBOT futures prices data from 1972 to 1997 confirms clearly that futures 
prices exhibit a seasonality. Another feature that is sugges t ed by the data is seasonal behaviour of the 
futures price volatilities. In this vein, Richter and SorensenI (2002) propose a model for the spot price of 
soybeans based on seasonal stochastic volatility. Geman and Nguven ( 2005l l also introduce a spot-based 
model for soybean prices with seasonality both for the price level and the (possibly stochastic) volatility 
level, back et ( 2013ll analyze data from corn, soybean, heating oil and natural gas markets and 
compare various spot-based models with deterministic seasonal volatility. They conclude that a volatility 
with s easonality is an important feature when valuing options on futures in these markets, back et al 


(2011) also study a futures-based model with seasonal stochastic volatility, which is essentially the 


Heston 


19931) stochastic volatility mod el with deterministic , seasonal mean-reversion level in the square-root 


process followed by the variance. Schmitz et al. ( 2013ll study calendar spread options in agricultural grain 
markets relying on a joint Heston model for the two underlying futures contracts. These two contracts 
share the same variance process, which has a const ant mean-r e versio n level, and therefore does not display 
seasonality. In the context of interest rate s , the Cox et al.l (198^ f CIR) model has been extended to 
time-dependent parameters by Maghsoodi ( 1996ll . and the HestonI (Il993ll model with time-dependent 
parameters, including the correlation between the spot price and its variance, has been studied by 
Ben^mou et al.l ( 201(TI1 . Let us also note that in the context of electricity markets, Lucia^^nd^,Schwar^ 
20021 1 give a detailed justification of the choice of seasonality function, as do iGeman and Roncoroni 


In parallel to his remark about seasonality in futures prices, Sorensen ( 2002ll confirms the Samuelson 
( 1965ll hypothesis that “the variations of distant maturity futures are lower than nearby futures prices.” 
We call thi s pattern the Samuelson effect. Popular futures-based models that incorporate this effect 
are those of Clewlo^jn^_StricMand (Il999al lbl) . The volatility function s used in these models are deter - 
ministic. Schneider and TavinI ( 201^ extend the multi-factor model of IClewlow and Strickland ( 1999b ) 
to incorporate stochastic volatility. Not only is stochastic volatility an incontestible empirical feature 
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of prices in futures markets, its inclusion also allows to calibrate the model to option volatility smiles 
and skews typically seen in futures option markets. In agricultural markets, a reflection of stochastic 
volatility is the recent introduction of several volatility indices on the CBOE/CBOT: the Corn Volatility 
Index (CIV) and Soybean Volatility Index (SIV) were introduced in 2011, and the Wheat Volatility Index 
(WIV) was introduced in 2012. 


In this paper, we extend the model introduced in Schneider and Tavinl ( 201511 to incorporate seasonal 
trends in the stochastic volatility processes. To achieve this, we begin by studying the mathematical 
conditions to impose on the seasonality function to guarantee that the generalized CIR process retains 
important features, such as existence and uniqueness of a strong solution, and positivity. It turns out 
that these conditions are very mild. These conditions appear to be not only interesting from a theoretical 
point of view, but also useful in practice, because different markets may need to be modelled with different 
seasonality patterns for the volatility. 


We then introduce the mod el with seasonal volatility and show how, by a generalisation of the results 


Schneider and Tavinl (|2015ri . the joint characteristic function of the log-returns of two futures prices 


can be obtained. It turns out that the Riccati ODE for the first function A is not affected, and only the 
integral ODE for t he second function B depen ds on 6 and is altered. Therefore, the same closed-form 
solution for A as in Schneider and Tavinl ( 2015h can be used. 

Next, we propose several specifications of seasonality functions and compare them. 


Then, we calculate implied volatility surfaces in our model. We also calculate calendar spread option 
prices and examine the effect of changes in the seasonality function on these prices. 

Regarding correlations, we study the effect seasonality has on the instantaneous correlation between 
two futures contracts in the case when the variances are deterministic. We find that the influence of the 
seasonality on the correlation depends on the magnitude of the exponential damping of the volatility 
factors. And from the calendar spread option prices we calculate implied correlations, where again we 
can observe how the seasonal pattern of the variance translates into a seasonal pattern of the implied 
correlation. 


The rest of the paper proceeds as follows. In Section[2]we discuss the CIR process with time-dependent 
drift. In Section[3]we define the proposed model and calculate the associated joint characteristic function. 
We also give the methods to be used for vanilla and calendar spread option pricing. In Section |4] we 
review various seasonality functions that can be used to specify the proposed model. Section [5] gives a 
numerical illustration of the implied volatility patterns produced by the proposed model with different 
seasonality functions. Section [5] deals with the seasonal behaviour of the instantaneous correlation when 
volatility is seasonal. In this section we also study the effect of seasonality on calendar spread option 
prices. Section [7] concludes. 


2 The CIR Process with Time-Dependent Drift 


To our knowledge, IHuII and Whitel ( 199f)ll were the first to consider extending the Cox et al. (1985) 
(CIR) interest rate model to time-dependent coefficients. They conclud e that in thi s gene ral case, it is 
no longer possible to obtain European bond option prices analytically. Maghsoodl (1996) also studies 
the “extended” CIR process in which the parameters k, ff and tr are time-dependent and finds, under 
certain conditions, the unique strong solution to the SDE describing the evolution of the process. 

In the context of the HestonI ( 199.811 stochastic vo latility model, the CIR p rocess represents the variance 
process of a stock price or foreign-exchange rate. Benhamou et al. ( 2010l) study the “time dependent 
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Heston model” and derive analytical formulas approximating European option prices. In their setup, the 
mean-reversion parameter k is constant, but the parameters 6, a and p (giving the correlation between 
the stock price, or foreign-exchange rate, and its variance) are all allowed to vary with time t. 

In the model introduced here, we only let the mean-reversion level given by 9 depend on time, while 
the other parameters k > 0 and cr > 0 (and later also p) remain constant. 

Let (f2,.A, P, be a filtered probability space, and let B = {Bt)t>o be a Brownian motion on this 
space. Let B = {ti,i = I,...} be a set of times having only finitely many points in every bounded interval, 
and let Z = {0 < < ^2 < •■• < < •■•} be the partition of M))" defined by T. Finally, let the seasonality 

function 9 : R))" R.’*' be piecewise continuous with respect to Z, and assume that it is bounded from 

below and above by positive constants 9 min and 9max- 

We will compare two processes v (seasonal) and v (non-seasonal), which are given, respectively, by 
the SDEs 

dv{t) = K {9{t) — v{t)) dt + a\lv{t)dB{t), (I) 

dv(t) = K {9min - v{t)) dt + ayjv{t)dB{t), (2) 

with identical parameters k > 0 , cr > 0 and initial conditions 0 < ■0(0) = vq < v(0) = vq. 

It is well known that ([2|) has a unique strong solution. The following result describes the solution to 

(mi- 

proposition 2.1 Assume that the seasonality function 9 is pieeewise continuous w.r.t. the partition Z 
and bounded by positive constants 9min o,nd 9max, *•£. for all t > 0,0 < 9min < 9{t) < 9max- Let 
the processes v and v be given by dH) and (0, respectively. Then: 

(i). The process H]) has a unique strong solution with continuous sample paths. 

(ii). P [Ot < ?;t, Vt > 0] = 1. 

(in). If the Feller condition < 2K9min is satisfied for 9min, then the process v is strictly positive. 

We prove this result in appendix 

Note that if the Feller condition is violated, then 0 can possibly reach 0, but it still cannot become 
negative. 

The piecewise continuity condition on 9 means that even discontinuous specifications of the mean- 
reversion level, such as the sawtooth function given by 

9{t) = a + b{t-to - [t-to\), (3) 

with a,b > 0 and to G [0,1[, pose no problems. Some other examples of the form of 9 are discussed in 
Sectional 

3 A Model with Seasonal Stochastic Volatility for Agricultural 
Futures 

3.1 The Financial Framework and the Model 

We begin by giving a mathematical description of our model under the risk-neutral measure Q. Let n > I 
be an integer, and let Bi,..., B 2 n be Brownian motions under Q. Let Tm be the maturity of a given 
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futures contract. The futures price F{t,Tjn) at time t,0 < t < Tm, is assumed to follow the stochastic 
differential equation (SDE) 

dF{t,Tm) = F(0,T™) = F^,o > 0. (4) 

i=i 

The processes Vj,j = are stochastic variance processes with time-dependent seasonal mean- 

reversion level assumed to follow the SDE 

dvj {t) = Kj {Oj (t) - Vj {t)) dt + <7j^Jvj{t)dB„+j (t), vj (0) = Vj^ > 0. (5) 

Various possibilities of the specification of the seasonal mean-reversion level functions 9j : Rj —>■ K+ 
are presented and discussed in Section 01 Note that the initial futures curve F{0,Tjn),m = 1,2,..., 
is exogenous in our model and can therefore accommodate any seasonal pattern shown by the futures 
prices. 

Eor the correlations, we assume 

{dBj{t),dBn+j{t)) = pjdt,-l < pj < l,j = l,...,n, (6) 

and that otherwise the Brownian motions Bj,Bk,k ^ j,j + n, are independent of each other. As we 
will see, this assumption has as a consequence that the characteristic function factors into n separate 
expectations. 

Eor fixed T^, the futures log-price In F{t,Tm) follows the SDE 

d\nF(t,T^) = , lnE(0,r„) = InF^.o- (7) 

Integrating o from time 0 up to a time T,T < Tm, gives 

n p'Y' - -| n n'Y' 

lnF(T,r„)-lnE(0,r„)= V / - - V / (8) 


We define the log-return between times 0 and T of a futures contract with maturity Tm as 

'F{T,T„ 


Xm{T) In 


F{Q,Tm) 


In the following, the joint characteristic function (j) of two log-returns Ai(T), X 2 {T) will play an important 
role. Eor u = (ui, U 2 ) € C^, (j) is given by 


4i{u) = (j){u\T,Ti,T2) =: 


exp 


The joint characteristic function $ of the futures log-prices lnE(r,ri),lnF(T,T 2 ) is then given by 

$(«) = exp ^i^UfclnE(0,rfc)^ •</'(«)• 


(9) 


( 10 ) 


Note that futures prices in our model are not mean-reverting, and that the log-price \nF(t^Tm) at time 
t and the log-return lnE'(T, T^) — in F{t,Tm,) are independent random variables. 

In the following proposition, we show how the joint characteristic function (j), and therefore also the 
single characteristic function c^i, is given by a system of two ordinary differential equations (ODE). 
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Proposition 3.1 The joint characteristic function (j at timeT < Ti,T 2 for the log-returns Xi{T), X 2 {T) 
of two futures contracts with maturities Ti,T 2 is given hy 


(j){u) = 4 )[u\T,Ti,T2) 

]^exp 0) (vj{0) + exp {Aj{0,T)vj{0) + Bj{0,T)). 


j=i 


where 


fjp{u,t) ='^Uke fj, 2 {u,t) = ^Uke 




/ \ ^7 '^7 

qj{u,t) = ipj - 




^j.T = 


Jo 


x^t 




and the functions Aj : {t,T) i— Aj{t,T) and Bj : {t,T) i— Bj(t,T) satisfy the two differential eguations 

dAj ^ 2 a2 

— KjAj + 2 '^j^j + dj = 0, 

- + Kj6j{t)Aj = 0 , 


dt 


with A,{T,T) = i^f^,,{u,T), Bj{T,T) = 0. 


The single eharacteristic function (f>\ at time T < Ti for the log-return Xi(T) of a futures contract 
with maturity Ti is given hy setting U 2 = 0 in the joint characteristic function. 


Note that the integrals Oj^r only depend on the specification of the seasonality fnnctions 9j and 
the maturity T. Therefore, their value can be calculated once and then stored, avoiding recalculations 
during repeated calls to the characteristic function. Of course, if 0 is a constant function, then the joint 
characteristic function given above is the same as the one given in Schneider and TavinI ( 2015 1. 

We prove this result in appendix 


3.2 Pricing Vanilla Options 


European opt ions on futures contracts can be priced using the Fourier inversion technique as described in 
Heston ( 1993 ) and Bakshi and Madan ( 2000ll . Let K denote the strike and T the maturity of a European 
call option on a futures contract F with maturity Tm > T. The function needed for this technique is the 
single characteristic function $1 of the futures log-price lnF(T, T^), given by $i(m) = 
with (fi (u) obtained from Proposition 13.11 European put options can be priced via put-call parity 
C-P = e-^^ iF{0,Ti)-K). 


3.3 Pricing Calendar Spread Options 

Calendar spread options (CSO) are popular options in commodity markets. To give a recent example, in 
February 2015 the Minneapolis Grain Exchange (MGX) introduced North American Hard Red Spring 
Wheat (HRSW) GSOs for trade on GME Globex. Calendar spread options are defined as follows. 
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Let two futures maturities Ti,T 2 , an option maturity T, and a strike K (which is allowed to be 
negative) be fixed. Then the payoffs of calendar spread call and put options, CSC and CSP, are 
respectively given by 

CSC{T) = {F{T, Ti) - F(T, T^) - K)+ , (11) 

CSP{T) = {K - {F{T, Ti) - F{T, T2)))+ . (12) 

To evaluate such options with a pricing model, the discounted expectation of the payoff must be calculated 
in the risk-neutral measure. As for “vanilla” European options, there is a model-independent put-call 
parity for calendar spread options: 


C'S'C(O) - CSPiO) = e-”-' (F(0, Ti) - F(0, T 2 ) - K). 


(13) 


Caldana and Fus^ ( 2013h show how to calculate CSC prices with models for which the joint characteristic 


function is known. Strictly speaking, their methods give a lower bound for the calendar spread option 
price, but usually this bound is very close to the true price. Note that in case the strike K = 0 the formula 
is exact. When we calculate CSO prices, we do so with this method . An alternati ve method, based on 
the 2- dimensional EFT algorithm, is given bv iHurd and Zhaiil ( 201(lll . We refer to Schneider and TavinI 
1)20151) and references therein for more details on calendar spread options. 


4 Seasonality Functions 

Below we present four types of seasonality functions that can be used as parametric forms to model 
seasonal variations of the volatility. For a given factor, 9 represents the seasonality term of the volatility 
dynamics and 9 is an integral involving 9 appearing in different expressions such as the characteristic 
function. For T > 0 and A € M, 0 is written 

9TiX) = [ 9{t)e^*dt. (14) 

Jo 

The presented seasonality functions 9 are parametric and work with three parameters: a, 6 and to- 
Parameter a controls for the volatility level. Parameter b controls the magnitude of the seasonality 
pattern and to corresponds to the time of the year when the volatility reaches its maximum. 

It seems relevant to consider different seasonality patterns as the reasons underpinning the seasonality 
phenomena in volatility may vary from a market to another. The first two patterns considered below are 
smooth and are based on the sinus function. The three others have points of non-differentiability and/or 
discontinuity, and may be used to represent a less regular evolution of the volatility. 

The sinusoidal pattern is written, with a, & > 0 and to € [0,1[, 

9{t) = a + bcos{2TT {t — to)), (15) 

„ bc^'^ 

= ,n , . n (27rsin(27r(r - to)) + A cos (27r(T - to))) 

A^ -I- 47r^ 

+ ^2 + (27rto) - A cos (27rto)) -f j (e^'^ - l) . (16) 

For a proof of this expression, we refer to Appendix 

The exponential-sinusoidal pattern is written, with a,b > 0 and to € [0,1[, 

0(t) = aexp (6cos (27r (t — to))). (17) 
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This parametric form for 9 is used in iBack et alJ (l201lli . There is no closed form expression for 0 t(A). 
The sawtooth pattern is written, with a,b > 0 and to € [0,1[, 


9{t) = a + b(t - to - [t - toj): 
^t(A) = y ^ - 


A 

5gAto 


^-AT 


a b \ T — — io 


[T-toj 


[T - toj ^ - 1 , 


(18) 


(19) 


k=l 


where [.J denotes the floor function, I is the indicator function and, by convention, X)fe=i = 0 if 
p < 1 . The proof leading to this expression is found in Appendix 1X1 


The triangle pattern is written, with a, 6 > 0 and to G [Oj l[i 


0 (t) = . 


- - (t-to - 


( 20 ) 


BriX) = j - 1 ) + 


5 gAto 

A 

2 a a 

-62 +226^ -21 1 ^ 


Z 2 +y-^e 2 +e “(^2 - to) j - e “( 2:2 - to)I{t^<i} 

2 


gAfc gAn 


-6 2 - 236^“ 
A 


{a>i}+^3e “I|„<i|-2i 


nT>to} 


+ (22 + r - to) - 22 ) 1 o[} 

(22 + T - to) + 22 'j I{j._t^g[_i _i [} 


-'r ^ 


( 21 ) 


with n = [T — toJ, a = r — to — [T — toJ, 2 i = A + 22 = ^ — y and 23 = 21 — a and with convention, 

X)fc=o = 0 if p < 0. The proof leading to this expression is found in Appendix lAl 


The spiked pattern is written, with a, & > 0 and to G [0,1[, 

2 


9{t) = a + b 


- 1 


( 22 ) 


_1 + |sin(7r(t - to))I 

This parametric form for 9 can be found in Geman and Boncoronil (|200fill where it is used to model the 
time varying intensity of a jump process. There is no closed form expression for 9t{X). 

Figure [1] presents the plots of these seasonal patterns with to = ^. 


5 Implied Volatility Smiles and Term-Structures 


In this section we present implied volatility smiles and term-structures produced by a one-factor version 
of the proposed model with seasonality. We have chosen the one-factor model here because our purpose 
is to illustrate the effect of seasonality on option prices. The considered options are vanilla options on 
futures. Following the market convention, their maturity is the same as the maturity of the underlying 
futures contract. The parameters of the considered model with seasonality are gathered in Table [TJ 
Only the parameters corresponding to the seasonality pattern change from one setting to another. The 




















seasonality functions 6 obtained with these parameters are those shown in Figure [TJ This numerical 
application is presented for an illustrative purpose and we take an initial futures curve that is flat in 
price at 100 USD. 


Table 1: Model parameters for different specifications of the seasonality functions used in the numerical 
illustrations. 


parameters 


seasonality function 


sinusoid 

exp-sinusoid 

sawtooth 

triangle 

spiked 


O.IO 

0.10 

0.10 

0.10 

0.10 

A 

1.00 

1.00 

1.00 

1.00 

1.00 


0.80 

0.80 

0.80 

0.80 

0.80 

a 

1.20 

1.20 

1.20 

1.20 

1.20 

p 

-0.25 

-0.25 

-0.25 

-0.25 

-0.25 

a 

0.25 

0.20 

0.10 

0.10 

0.10 

b 

0.15 

0.68 

0.30 

0.60 

0.30 

to 

7/12 

7/12 

7/12 

7/12 

7/12 


In Figures [2] and [3] we present, for the different seasonality patterns, the implied volatility smiles 
obtained for different maturities and the term-structure of implied volatility for at-the-money options. 
The obtained term-structures exhibit both seasonality and the Samuelson effect. The obtained strike- 
structures exhibit smiled shapes. The implied volatility produced with the sinusoidal and exponential- 
sinusoidal patterns are similar to each other. For the sawtooth, triangle and spiked patterns, the irreg¬ 
ularities of the function 0 seem not to be transferred to the implied volatility. It can also be observed 
that the choice of the seasonality pattern seems to have very little impact on the shape of the volatility 
smile. This last remark is particularly striking if one compares the volatility term-structures produced 
by the sinusoidal and triangle patterns. 


6 Seasonal Stochastic Correlation and Calendar Spread Option 
Prices in the Multi-Fact or Model 

6.1 Seasonal Stochastic Correlation 

We will show in this section that if we specify our model with two or more volatility factors, then the 

correlation of the returns of two given futures contracts is influenced by the seasonality functions. In 

other words, the correlation becomes seasonal. 

Let us take the 2-factor model for an illustration. Futures returns follow the SDE 

^ \^t -^m) 

and the two variance processes follow the SDEs 

dvi{t) = Ki {9i{t) - vi(t)) dt + <TiyJvi{t)dB 3 (t), (24) 

dv2{t) = K2 {d2{t) — V2{t)) dt + cr2^Jv2{t)dBi{t). (25) 

The correlations are given by {dBi{t), dB^it)) = pidt, {dB 2 {t),dBi{t)) = p 2 dt, and all other correlations 
are zero. 
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Define 


(26) 


V^j{t) := ( 


dF{t,T,) 

F{t,Ti) 


dF{t,T,) 

F{t,T,) 


)/dt. 


Then the instantaneous correlation p{t) at time t is given by: 


P{i) 


Vl2{t) 

\/Vii{t) y^V22{t) 


(27) 


Using (1^ together with (1^ gives for the instantaneous correlation 


P{t) 


e A2(Ti+T2 

^g-2Ai(Ti-t)yi(^) _|_ e-2A2(ri-i)g2(t)^e“^'^i(^2-t)^j^(i) _|_ g-2A2(T2-i)g2(t) 


(28) 


In contrast to the 1-factor model, the instantaneous correlation p(t) in the n-factor model, with n > 2, 
is stochastic. 

To illustrate the seasonality of the correlation function p, we consider the case cti = (72 = 0 when 
both variances vi and V 2 are deterministic functions of time t. The two SDEs ([S]) then become ordinary 
differential equations 

v'jit) = Kj {0j(t) - Vj{t)) , Vj{0) = Vjp > 0. (29) 

It is straightforward to see that the solution to (l29)) is given by 


Vj{t) 


= e ^ 



= e ^ 


.Vjp + Kj§t{K^ 


(30) 

(31) 


Note that the same transform function 9t{K) = J* e'^^9{s)ds already introduced in (1141) appears again. 
For several specifications of 0, the function 9t is available in closed form, and it is easy to plot the 
correlation function p given by ([28]). 

Figure m presents the instantaneous correlation obtained with the two-factor model with seasonality 
in the special case cti = tT 2 = 0. It also plots the correlation obtained with the version without seasonality 
and the other terms involved in expression (E51) . For the seasonal model, only the first factor is seasonal 
and follows the sinusoidal pattern m- To produce Figure 01 we consider two cases in order to illustrate 
different seasonal patterns of the correlation when compared to the non-seasonal case. The non-seasonal 
case is our benchmark, which is obtained by setting & = 0 in the seasonality function of each factor. The 
model parameters for these cases are gathered in Table 0J 

It is worth remarking that, compared to the non-seasonal case, the instantaneous correlation is 
modified by the presence of seasonality in the dynamics of the first factor. In the first case (high A for 
the seasonal factor) this correlation is lower at the beginning of the period under scrutiny and higher at 
the end. In the second case (low A for the seasonal factor) it is the opposite. The reasons underpinning 
the influence of Ai and A 2 on the seasonal behaviour of the correlation are still conjectural and need to 
be further investigated. So far, it seems on the one hand that if we have A 2 < I < Ai, then the effect 
of the first seasonal function on the instantaneous correlation is in the opposite direction, i.e. during 
periods of high variance the correlation is decreased, and during periods of low variance the correlation 
is increased. This pattern can be observed in the three l.h.s. panels of Figure 0] illustrating Case 1. On 
the other hand, if we have Ai < I < A 2 , then the effect of the first seasonal function on the instantaneous 
correlation is in the same direction, i.e. during periods of high variance the correlation is increased, and 
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Table 2: Model parameters in the two cases used to illustrate the seasonal behavior of the instantaneous 
correlation in the two-factor version of the model. 


parameters 

Case 1 

Case 2 

I’Ol 

0.10 

0.06 

f02 

0.04 

0.04 

Ai 

2.00 

0.50 

A 2 

0.50 

2.00 

Kl 

1.00 

1.00 

K2 

1.00 

1.00 


0.00 

0.00 

^2 

0.00 

0.00 

Pi 

0.00 

0.00 

P2 

0.00 

0.00 

Oi 

0.10 

0.06 

02 

0.04 

0.04 

hi 

0.09 

0.05 

b2 

0.00 

0.00 

toi 

0.00 

0.00 

t02 

0.00 

0.00 


during periods of low variance the correlation is decreased. This pattern can be observed in the three 
r.h.s. panels of Figure S] illustrating Case 2. 

Note also that the instantaneous correlation is not always decreasing with t. This remark holds true 
for the model with seasonality as well as for the case without seasonality. Figure [5] presents a non¬ 
decreasing instantaneous correlation obtained with the two-factor model with and without seasonality in 
the special case cti = tT 2 = 0. The other parameters are chosen such that the curves are non-decreasing. 
In the presented illustration the correlation seems to be constant at first sight, but on closer inspection 
can be seen to be decreasing then increasing. 


6.2 The Effect of Seasonality on Calendar Spread Option Prices 


In this section, we investigate the effect of seasonality on the prices of calendar spread options. We do 
so by means of a numerical study in which we consider different cases of model parameters and compute 
calendar spread option prices with and without seasonality. 


This investigation is conducted with the two-factor model with stochastic volatility on both factors 
but seasonality only on the first one in order to isolate the effect of seasonality. The form of the seasonal 
pattern i s again the sinusoidal one given by equation (fTbll . Calendar spread option prices are obtained 
with the ICaldana and Fnsail (2013) method that is well suited for our model. Implied correlations are 
extracted using a numerical root search. 


Table [3] presents, for different magnitudes of seasonality, calendar spread option prices obtained with 
the two-factor model. Case 1 corresponds to the model without seasonality, Case 2 corresponds to a 
moderate seasonality and Case 3 to a stronger seasonality. The model parameters for these three cases 
are shown in Table m Figure |6] plots the implied correlation term-structures obtained from calendar 
spread options in the considered cases. 
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Prices presented in Table [3] show the effect of seasonality on calendar spread options. Irrespective of 
the strike, prices are increasing with the seasonality magnitude for the two maturities before and at a 
multiple of to and decreasing for the two maturities after a multiple of to- 

Figure [5] shows that seasonality has a noticeable effect on the implied correlation term-structure 
for calendar spread options. Compared to the non-seasonal case, implied correlations are higher for 
maturities before a multiple of to and lower for maturities after a multiple of to- 


Table 3: Calendar spread option prices obtained with the two-factor model with stochastic volatility and 
different magnitudes of seasonality on the first factor. The reported prices are for different strikes and 
maturities. The difference between the underlying futures maturities is held constant at 6 months. The 
corresponding model parameters are shown in Table SI 


T 

T1 

T2 

Case 1: 

no seasonality 

Case 2: moderate seasonality 

Case 3: 

strong seasonality 

K = -10 

o 

II 

K = 10 

K = -10 

o 

II 

o 

II 

K = -10 

o 

II 

K = 10 

0.33 

0.33 

0.83 

10.4381 

2.5993 

0.4951 

10.5605 

3.0931 

0.6822 

10.6647 

3.4190 

0.8351 

0.58 

0.58 

1.08 

10.6718 

3.4672 

0.9188 

10.8369 

3.7314 

1.1026 

10.9583 

3.9252 

1.2400 

0.83 

0.83 

1.33 

11.1648 

4.5229 

1.6366 

11.1425 

4.3690 

1.5636 

11.1291 

4.2429 

1.5120 

1.08 

1.08 

1.58 

11.3550 

4.6534 

1.8404 

11.2930 

4.5533 

1.7709 

11.2459 

4.4773 

1.7187 

1.33 

1.33 

1.83 

11.1894 

4.1480 

1.6234 

11.2674 

4.4534 

1.7695 

11.3323 

4.6612 

1.8843 

1.58 

1.58 

2.08 

11.1872 

4.4875 

1.7783 

11.3086 

4.6726 

1.9165 

11.3983 

4.8091 

2.0197 

1.83 

1.83 

2.33 

11.5326 

5.2374 

2.3126 

11.4832 

5.0666 

2.2079 

11.4488 

4.9289 

2.1312 

2.08 

2.08 

2.58 

11.6266 

5.2422 

2.3840 

11.5474 

5.1193 

2.2922 

11.4873 

5.0262 

2.2233 

2.33 

2.33 

2.83 

11.3966 

4.6788 

2.0671 

11.4623 

4.9321 

2.1982 

11.5160 

5.1052 

2.2992 

2.58 

2.58 

3.08 

11.3529 

4.8929 

2.1509 

11.4589 

5.0521 

2.2717 

11.5367 

5.1692 

2.3617 

2.83 

2.83 

3.33 

11.6528 

5.5345 

2.6146 

11.5933 

5.3604 

2.4988 

11.5511 

5.2210 

2.4132 


7 Conclusion 

We introduce a multi-factor seasonal stochastic volatility model for futures contracts that is capable of 
reproducing the Samuelson effect. We show that the model can accommodate very general specifications 
of the seasonality functions, including only piece-wise continuous ones. As an illustration, we suggest 
five different seasonality functions, some of which are familiar from the literature, and provide details 
of how to incorporate these into the joint characteristic function of the model in a numerically fast and 
efficient way. 

In a series of examples, we show that this model can reproduce seasonal implied volatility surfaces 
of European options on futures contracts. Furthermore, implied correlations calculated from calendar 
spread option prices also show seasonal patterns. Finally, we demonstrate that the instantaneous correla¬ 
tion between the returns of two futures contracts with different maturities is both stochastic and seasonal, 
and make a conjecture about the relationship between the magnitudes of the Samuelson damping factors 
and the effect of the seasonality functions on this correlation. 
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Table 4: Model parameters in the three cases used to compute calendar spread option prices and the 
corresponding implied correlations. 


parameters 

Case 1 

Case 2 

Case 3 

^^01 

0.10 

0.10 

0.10 

V02 

0.04 

0.04 

0.04 

Ai 

2.00 

2.00 

2.00 

A2 

0.50 

0.50 

0.50 

Ki 

0.80 

0.80 

0.80 

K2 

0.80 

0.80 

0.80 

CTl 

1.20 

1.20 

1.20 

02 

0.90 

0.90 

0.90 

Pi 

-0.25 

-0.25 

-0.25 

P2 

-0.25 

-0.25 

-0.25 

«! 

0.25 

0.25 

0.25 

02 

0.10 

0.10 

0.10 

hi 

0.00 

0.15 

0.35 

^2 

0.00 

0.00 

0.00 

toi 

7/12 

7/12 

7/12 

to2 

7/12 

7/12 

7/12 


A Proofs 


In this appendix we prove Propositions 12. II and 13. II We also give the proofs of expressions (ITB)) and dTiH) . 
Proof of Proposition [27ll 


(i). The drift function b(t,Vt) '■= K{9{t) — v{t)) in ([T]) is Lipschitz continuous w.r.t the second argument, 


i.e. 


\b{t,x) -b{t,y)\ < K\x-y\, 


where we can choose K = k, since \b(t, x) — b{t, y)\ = k\x — y\. Then Proposition 2.13 (Yamada and 
Watanabe) of Karatzas and Shrevd (|l988ll guarantees the existence of a unique strong solution to 
o with continuous sample paths. 


(ii). 


The comparison result given in Proposition 2.18 of Karatzas and Shrevel (1988) establishes Vt > Vt 
a.s. for all t > 0 under the hypothesis that the drift function b{t,vt) is continuous. Now, if 9 
has a discontinuity at time ti, we know from this argument applied to the interval [0,ti[ that 
Vt < ftVt € [0,ti[ (a.s.). It then follows from the continuity of the sample paths that 
(a.s.), and we can apply the argument again to the interval ]ti,t 2 [ to obtain vt < VtVt G]ti,t2[ 
(a.s.). Since by assumption the set T of times where 9 has discontinuities has no limit points, we 
can proceed in this manner to cover all of. 


1 + 


(hi). The Feller condition cr^ < 2K9min for 0mm implies the strict positivity a.s. of v. The strict positivity 
of V itself therefore follows immediately from (ii). 


□ 


Proof of Proposition r3.ll The proof is an extension of the proof of Proposition 2.1 of ISchneider and Tavin 


(|2015ll to the case where the variance mean-reversion level 9 is time-dependent. Going from 9 to 9{t) 
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leads to changes in two places. The first is in Lemma A.l of ISchneider and TavinI (|2015r) . which needs 
to be modified as follows. 

Lemma A.l Let 9 : ^ K’*' be the seasonal mean-reversion level function, and let 


0j,{X) := / e^*9{t)dt 

Jo 


be its transform. Then 

rT 


10 


fi(t)y/v(t)dB{t) = [fi(t)v{t)]Q - /i(O)k0t(A) + (k - A) / fi(t)v{t)dt. 


(32) 


Proof. Multiplying equation ([5]) by fi{t) and then integrating from 0 to T gives 


fi{t)dv{t) = [ fi{t)K{eit) - v{t))dt + a [ fi{t)y/v(fydB{t). 
Jo Jo 


Using Ito-integration by parts Isee l0ksendall (2003)), we also have 

[ fi{t)dv{t) = [fi{t)v{t)]g - [ v{t)dfilt) 

Jo Jo 

= [hitMt)]o - ^ i fi{t)v{t)dt. 

Jo 

Equating the right hand sides of equations (IM)) and (l34ll gives 

O' / h(.t)\/W)d'B{t) = [fi{t)v{t)]^ - X [ fi{t)v{t)dt- [ fi{t)K{e{t) - v{t))dt 
Jo Jo Jo 

pT rT 

= - K fi{t)e{t)dt-\-{K-X) fi{t)v{t)dt 

Jo Jo 

.T pT 

= - fiiO)K + (k - A) / fi{t)v{t)dt 

Jo Jo 

= [/i(0'o(0]o - /i(O)k0t(A) + (k - A) / fi{t)v{t)dt, 


(33) 


(34) 


which proves the lemma. 


□ 


The second change i n the proof is due to the appearance of 9 in the generator of the process v. As 


Schneider and Tavin ( 2015h . let the function h be given by 


h(t, u) = E 


exp (i^fi{T)v{T) + q{s)v{s)ds 


Now h satisfies the PDE 


— {t,v)+K{9{t) 


r)h 1 h 

'o(0)-5-(i)^) + + <l{t)v{t)h{t,v) = 0, 


dv 


dv^ 


(35) 


with terminal condition 

h{T,v) = exp (^i^fi(T)v{T)^ . 
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( 36 ) 


Again, we know from DufRe et alJ (200 Clli that h has affine form 

h(t, v) = exp {A(t, T)v{t) + B(t, T )), 
with A(T,T) = i-^fi{T),B(T,T) = 0. Putting (1361) in (1351) gives 

Bt + Atv + K{9{t) — v)A + ^a^vA^ + qv = 0, 
and collecting the terms with and without v leads to the two ODEs 

At — kA + —^ = 0, 


Bt + Ke{t)A = 0. 


This completes the proof of the proposition. 


(37) 

(38) 

□ 


Note that 6 only appears in the second ODE 


PRI) , and that therefore the closed-form expression 
previously given for A in Schneider and TavinI ( 2015ll can still be used. Only the function B changes due 
to the time-dependence of 6. 

Proof of expression (I16|) . Transform of the sinusoidal pattern. 

With y = t — to 


0t(A) = / {a + b cos {2Tr {t — to))) e^^dt 

Jo 

= Y ~ l) + be^*° [ cos {2'!Ty)e^^dy. 
A J-tn 


A primitive of ?/1—>■ cos {2TTy)e^^ is 




y 


A^ -I- dTT^ 


(A cos {2TTy) + 27r sin (27rj/)), 


(39) 

(40) 

(41) 


and 



cos {2TTy)e^'^dy 


Av 


A^ -I- 471^ 
gA(T-to) 


(A cos (27ry) -|- 27rsin (27r?/)) 


T—to 


J —to 


A^ -I- 47r2 

g-Ato 

A^ -I- 4:TT^ 


(27rsin {2tt{T — to)) + A cos (27r(T — to))) 
(27r sin (27rto) — A cos (27rfo)) 


Proof of expression (I19|) . Transform of the sawtooth pattern. 

With T > 0 and to € [0,1[ 


^t(A) = / {a + b{t-to-[t-to\))e dt 


/ [a + b{t — to)) e^^dt — b / [t — to\e^^dt. 

Jo Jo 


(42) 

(43) 

(44) 

□ 


(45) 

(46) 
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The first integral is computed as 

rT 


j (a + 6 (i - io)) t Q ^ + 6 - t - ip 

The integral involving the floor function can be split, with y = t — to-, as 

f [t — to\e^*dt= f [y\ + f [y\ dy. 

Jo Jo J —in 


(47) 


(48) 


Noting that [yj = — 1 for y e [—to, 0[ we have 

fO 

'-to 

The other part of the term with the floor function can be written, when T > to, as 


r LyJ e^^y+^°^dy = y (1 - 6^*“) . 

J-tn ^ 


pT—to /n — 1 pk-\-l pn-\-a. 

/ [yJ e^^^^*°'^dy = I ^ k / e^^dy + n e^^dy 

J 0 Vfe 0*^^ 


D-^^0 


A 


^gA(n+a) _ gAfc 


gA(n+0!) _ gAn 


E< 

fe=i 


with n = [T — toj and a = T — to — [T — toj • 

When 0 < T < to, as to < 1, [T — foj = —1 and [y\ = —1 for y G [T — to, 0[ so that we have 

rT-to pXto 

[?/J e^^'^~^*°^dy = - 
lo 


(•T-to pXto , 


Gathering the components, the integral involving the floor function can now be written as 

J [t - toj [ [T - toj ‘o) _ / ^ g^fc j _|_ _l_ g ->'‘0 _ ^ 


(49) 


(50) 


(51) 


(52) 


Proof of expression (I21|l . Transform of the triangle pattern. 

With T > 0 and to G [0,1[ 


0 r(A)= j [a + b 


--{t-to- L^-M) 




= ^(e--l)+5 


--{t-to- [t-to\) 


e^*dt. 


□ 


(53) 

(54) 
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With y = t — to, the last integral becomes 


-{t-to- U-ioJ) 


e^*dt = 


' —to 


f^T—to 


-{y- bj) 2 “ “ LyJ) e^^dyj . 

(55) 

Two integrals remain to be computed, one on [—to,0] and the other on [0,T — to]. To compute the 
first, one needs to distinguish two cases. When to S [0, ^], it can be computed as 


l-to 


2 - b - bJ) 


e^^dy =j[z2- {Z2 - to) e 


and when to G]^, 1[, it is 


-*o 


-{y- bJ) 


e^ydy = ^ (^2 + + {z2 - to) 


(56) 


(57) 


1 _ 1 

2 A 

First, when T > we have 


with Z 2 = ^ — To compute the other integral, on [0,T — to], one needs first to distinguish two cases. 


n'T'—to 


- b - bJ) 


z^^dy = ^ 


1-1 


' 


- b - bJ) 


e^^dy - 


nn-\-<y. 


- b - bJ) 


e^ydy, (58) 


with n = [T — toj and a = T — to — [T — to] ■ For k = l,...,n — 1, the integral in the sum can be 
computed as 


fk+l 


- b - bJ) 


z^^dy = 


ffc+i 


- b - bJ) 




rk+l 


- b - bJ) 


f>k+ 


2 /I 


--y + k] e^^dy - 


/•Ai+l / I 

I ^ ( n - 2 / + ) e.^'^dy 






where -2^2 = 5 + a- becomes 


’^-1 rk+l 


E 

fc=0 ’ 


- b - bJ) 


e^ydy =j(je^+ Z2e^ “ E * 

^ ^ k^O 


.Afc 


The integral on [n,n + a] is computed, as 


^n+a 
J n 


\-{y-[y\) 


^An 


= — 


2:36 


Aa-ri 


><n^U 


+ ( -e= -zoe^°‘ ) -^1 




(59) 

(60) 
(61) 

(62) 

(63) 


with zo = zi — a. 

The second case is when T G [0,to[. In this case, the integral on [0,T — to] becomes 


nT — to 


^ - b - bJ) 


.^Vdy = i ((e^(^-*o) {z 2+T- to) - Z2) I{T-toe[- 1 . 0 [} 

~ ( A®~" {z 2 + T - to) + 2:2^ I{T-toG[-l -i[} j ■ 


Gathering the components now gives the result. 


( 64 ) 

□ 
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Figure 1: Upper left: sinusoidal pattern. Upper right: exponential-sinusoidal pattern. Center left: 
sawtooth pattern. Center right: triangle pattern. Lower: spiked pattern. 
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Figure 2: Upper panel, implied volatility smiles (left) and term-structure (right) obtained using the 
sinusoidal seasonality pattern with parameters in Table [1] Lower panel, implied volatility smiles (left) 
and term-structure (right) obtained using the exponential-sinusoidal pattern with parameters in Table 

m 


21 



























Figure 3: Upper panel, implied volatility smiles (left) and term-structure (right) obtained using the 
sawtooth seasonality pattern with parameters in Table [1] Center panel, implied volatility smiles (left) 
and term-structure (right) obtained using the triangle seasonality pattern with parameters in Table 
[TJ Lower panel, implied volatility smiles (left) and term-structure (right) obtained using the spiked 
seasonality pattern with parameters in Table [TJ 
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Figure 4: Upper panel: instantaneous correlation obtained with the two-factor model with seasonality 
{red) and without seasonality {blue) for the two cases detailed in Tabled (case 1, left and case 2, right). 
Center panel: terms involved in expression (12811 . 0i{t) in red, vi{t) in black and 02(i) in green, V 2 {t) in 
blue, for the two cases detailed in Table [2] (case 1, left and case 2, right). O2 and V2 are identical as 
the second factor does not have seasonality. Lower panel: difference between instantaneous correlations 
obtained with and without seasonality for the two cases detailed in Tabled (case 1, left and case 2, 
right). 
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Figure 5: Instantaneous correlation obtained with the two-factor model with seasonality {red) and with¬ 
out seasonality (blue). Parameters are chosen such that both curves are non-decreasing. 
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Figure 6: Term-structure of implied correlation from calendar spread option prices obtained with the 
two-factor model with stochastic volatility and different magnitudes of seasonality on the first factor. 
The different magnitudes considered are: no seasonality {black), moderate seasonality {blue) and strong 
seasonality {red). The corresponding model parameters are found in Table 01 
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